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A general approach for investigating transport phenomena in porous media is presented. This 
approach has the capacity to represent various kinds of irregularity in porous media without the need 
for excessive detail or computational effort. The overall method combines a generalized Effective 
Medium Approximation (EMA) with a macroscopic continuum model in order to derive a transport 
equation with explicit analytical expressions for the transport coefficients. The proposed form of the 
EMA is an anisotropic and heterogeneous extension of Kirkpatrick's EMA [Rev. Mod. Phys. 45, 574 
(1973)] which allows the overall model to account for microscopic alterations in connectivity (with the 
locations of the pores and the orientation and length of the throat) as well as macroscopic variations 
in transport properties. A comparison to numerical results for randomly generated networks with 
different properties is given, indicating the potential for this methodology to handle cases that would 
pose significant difficulties to many other analytical models. 
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I. INTRODUCTION 

There are a large number of applications for both nat- 
ural and designed porous materials including those in 
fuel cell technology physiological transport phenom- 
ena [5] , and heat exchangers [3J 0] . By natural porous 
media, we mean structures with random features; like ir- 
regular pores in view of their sizes and shapes as well 
as their interconnections. On the other extreme, the 
term designed porous media [5] refers to purposefully 
constructed porous structures where pores are designed 
to perform functions under certain constraints; with heat 
exchanger tube bundles in cross flow as a very good ex- 
ample |()1[7] . While porosity can be very high for designed 
porous materials, natural porous materials typically have 
much lower values for porosity [5]. 

Although it is typically the macroscopic transport 
properties of either a designed or a natural porous 
medium in which we are interested, these properties are 
irrevocably tied to the structure of the medium on a mi- 
croscopic scale. However, due to the complexity of the 
pore-scale geometry, full analytical results at a micro- 
scopic scale are nearly impossible, while numerical simu- 
lation tends to be highly computationally intensive and 
require excessive detail. The approaches of volume av- 
eraging [8] and conditional averaging [9] provide simpli- 
fied analytical equations, but these equations remain dif- 
ficult to associate with the microstructure of the porous 
medium, especially its connectivity. 

A convenient approach that allows at least some of 
these microscopic properties to be accounted for is pore 
network modeling, which represents the pore space as 
a large number of pores connected together by narrow 
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throats with simplified geometry for both. With the use 
of numerical simulations on networks reconstructed from 
real-world media, realistic predictive estimates of trans- 
port coefficients have been made [TD1 [TT] . However, such 
reconstruction processes are relatively difficult, and net- 
work simulations can be quite computationally intensive 
for very large numbers of pores. Analytical models with 
suitable averaging remain highly desirable for the pur- 
poses of efficiency and robustness. 

Unfortunately, there are relatively few analytical ap- 
proaches to modeling of porous media that take the net- 
work microstructure into account. One such approach is 
the continuous-time random walk (CTRW) method [T2T 
114] ; however, apart from some recent efforts [12] , the cur- 
rent CTRW formulation usually requires its key param- 
eter to be fitted to experimental data. Another group 
of approaches is the effective medium approximations 
(EMAs), which can be applied to stochastic networks, 
replacing them with deterministic ones on the condition 
that the average of the local fluctuations must be zero. 
Typically, an effective medium approximation is applied 
to a regular network structure [15] . and then transport 
coefficients are calculated by assuming the existence of 
a smooth macroscopic field, and calculating local fluxes 
and hence transport coefficients based on this Smooth 
Field Assumption (SFA) [TH1 E]. In general, effective 
medium approximations are reasonably accurate as long 
as the medium is not near its percolation limit |18j . 

The limitations of standard EMAs are demon- 
strated by network structures reconstructed from three- 
dimensional (3D) imaging of samples of real- world porous 
media [10 . It is clear that regular lattices are not rep- 
resentative of the complex, irregular nature of porous 
media, and although they can often be tuned in order 
to match empirical data, the parameters in this tuning 
are typically just fudge factors. In real-world media, the 
pores and throats vary at macroscopic scales, the throat 
properties are often anisotropic, and the coordination 
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number exhibits significant variance. Some extensions 
of the effective medium theory have been proposed that 
deal with one or two of these issues (anisotropic EMA 
[TM2"T] . effective throat area [55], variable coordination 
number |23p. but an effective medium theory which deals 
with truly irregular networks of the kind we expect from 
natural porous media seems to be absent from the liter- 
ature. 

In this paper, we propose a generalized single-bond 
EMA which does not require a recurring network topol- 
ogy, and can have both anisotropy and large-scale varia- 
tion. This is combined with an explicit derivation of the 
transport equations for the effective medium network, in 
the form of a diffusion-like partial differential equation. 

Section[n]describes a general pore-scale model of trans- 
port within porous media, in which a porous medium 
is represented as an irregular pore-throat network which 
does not conform to any lattice pattern. Pore-scale trans- 
port is then modeled using general transport equations 
based on conservation principles. 



Section III presents a generalized EMA, which is used 
to account for microscopic variations due to poor connec- 
tivity within the medium. The result of this approxima- 
tion is a discrete pore network, but one with deterministic 
throat conductances. 

Section |IV| describes the continuum approximation 
used in order to derive an analytical expression for the 
transport coefficients in the medium. 

SectionfVlshows a comparison of the analytical method 



of Sections |III| and IV against numerical simulations for 
randomly generated pore networks. 



II. PORE-SCALE MODEL 

In order to achieve a tractable representation of a 
porous medium at the microscopic level, we begin our 
analysis with a stochastic pore-throat network model, 
representing the medium as a network of pores connected 
together by narrow throats. Although the network model 
we consider is nonetheless a simplification of the struc- 
ture of real porous media, we avoid many shortcomings of 
network models by considering a general class of stochas- 
tic networks. Using a stochastic model allows for realistic 
modeling of porous media by only using the important 
data about the statistical properties of a medium at a 
microscopic level, such as types and distributions of con- 
nections and pores. The stochastic properties then serve 
to represent uncertainty about the exact microstructure 
while taking such information into account. Moreover, 
the types of networks considered are capable of repre- 
senting porous media with various kinds of irregularity, 
including macroscopic heterogeneity, significant variation 
in co-ordination numbers, and anisotropy. 

Based on a simplification of real media, we consider the 
pore throats to have negligible storage capacity (typically 
volume), so that accumulation of the transported quan- 
tity occurs within the pores. However, the throats are 



the primary consideration in determining the throughput 
within the medium, because overall transport is limited 
by the constricted nature of the throats. In particular, 
if there were no throats at all, the medium would con- 
sist only of isolated pores, and no transport would take 
place. As such, transport is primarily modeled within 
the throats, and at the pore/throat interfaces. 

We represent the medium by a network of N pores, 
which we enumerate with integers 1, . . . ,N, where N 
must be large in order for a distinction to emerge between 
macroscopic and microscopic properties of the medium. 
A single pore within this network, labeled by its num- 
ber z, is then described by its location in n-dimensional 
space, Xi = (xj, . . . , a;™), which we take to be the posi- 
tion of its center, and a vector of additional properties £i 
(e.g. pore volume). For our purposes, we combine the 
physical coordinates and the additional properties into 
the extended vector yi — taking values in some 

vector space S, which we will call "extended space". 

In order to describe the statistical uncertainty that is 
inherent in modeling the complex and disordered struc- 
ture of a porous medium, the extended coordinate vectors 
are considered to be random variables Yi, Yz, . . . , Yn, 
corresponding to the positions and properties of pores 
sampled from the medium without replacement. As such, 
any individual Yi has an identical marginal PDF (Prob- 
ability Density Function), fyiy). For notational conve- 
nience, whenever the random variables considered are ev- 
ident from the arguments, the subscripts will be dropped 

- e.g. f(y) for the single-pore distribution. 

Throughout this network, there are throats connect- 
ing pairs of pores - these are thin links through which 
transport between pores can occur. In general, like the 
pores, these throats could have various properties associ- 
ated with them. However, for the purposes of the pore- 
throat network model, these properties are combined into 
a single value which we call the conductance of the throat, 
by analogy to electric circuits. For two pores i and j the 
conductance of a throat connecting them is the value gij, 
which describes the ease of transport through the throat 

- a value of would mean there is no direct link, while 
a value of 00 would effectively combine the two pores 
into a single pore. Moreover, we do not consider the 
conductance of an individual throat to depend upon the 
direction of transport, i.e. — gu. 

As with the pores, we consider this conductance to be 
a stochastic variable, , since individual throats in the 
network could have a variety of values for conductance. 
It is reasonable to assume that this conductance depends 
primarily on the properties of the throat itself, and on 
the pore-throat interfaces on either side. Of course, we 
can expect correlation between the throat and pore prop- 
erties - for example, larger pores would generally have 
larger throats. With this in mind, we express the distri- 
bution for conductances for any given throat by 

= f{g\vuVj), (1) 
which accounts for the effect of pore properties and lo- 
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cations on the throat conductance, as well as for inde- 
pendent randomness within the throats themselves. It is 
important to note that throat length is also accounted 
for in this representation, since the physical coordinates 
of the pores are included in and y 3 . 

In a realistic porous medium, it is clear that the vast 
majority of pairs of pores do not have any throat con- 
necting them at all, especially when the pores are spa- 
tially distant. Considering this, we let ^(yi,yj) be the 
probability that a throat between pores with extended 
coordinates Vi,yj is present, and so the conductance 
distribution for the throat can be written as 

fig I vuVj) = ^{yr,yj)f*{9 I vuVj) 

+[l-*(tfc,|fc)]*(ff), (2) 

where /* is the distribution over nonzero values of con- 
ductance (i.e. conditional upon G > 0), and S is the 
Dirac delta function. 

This stochastic network model is a very general type of 
random graph, with the addition of random edge weights 
- the conductance values representing the transport prop- 
erties of the throats. Of the many classes of random 
graph studied in the literature, the closest well-known 
one is the Random Geometric Graph (RGG) [23] due to 
the spatial nature of the node parameters and localized 
connectivity. Another approach is htness-based connec- 
tivity |2S], but this concept is targeted towards the gen- 
eration of scale-free networks and hence is not used here. 

Having defined the throat conductance, we can derive 
a general transport equation valid for a variety of dif- 
ferent transport phenomena, including diffusion, conduc- 
tion heat transfer, flow of electrical current, transport 
across biological membranes, and linear cases of fluid 
flow (e.g. incompressible Newtonian Stokes flow). We 
consider inter-pore transport to be driven by a "poten- 
tial difference" between two pores, while the transport 
phenomenon represents the tendency for these potentials 
to equilibrate. 

In this paper, we will use the terminology of molecular 
diffusion, and so we identify the potential of a pore as 
the concentration of the substance within the pore, Cj = 
Trii/Vi, where rrii is the amount of the substance in the 
pore, and Vi is the volume of the pore. There is a net 
flow of particles from pores with higher concentrations to 
those with lower concentrations, resulting in a flow rate 
through any given throat (from pore i to pore j) defined 
by 

Qij = -Qji = 9ijiCi ~ Cj)- (3) 

With the use of a conservation law within the i-th pore, 
we can finally derive the master equation, 

dm N N 

= = g v ( 4 ) 

3=1 3=1 

which is a system of N differential equations that de- 
scribes transport within the medium at the pore scale. 



As noted previously, this equation is applicable to a wide 
range of transport processes; expressions for throat con- 
ductance g in terms of the geometric properties of the 
throat and the properties of the materials or substances 
involved are given in Appendix [X] for several of these 
transport phenomena. 

III. NONUNIFORM EFFECTIVE MEDIUM 
APPROXIMATION 

If the porous medium is very well-connected, global 
transport can be directly approximated using the contin- 
uum approximation described in Section |IV| using the 
mean or expected value for the conductance of each 
throat, 

poo 

g{yi,yj)= gfig\vhVj)dg. (5) 

Jo 

However, in general, the media we consider do not meet 
this requirement. In order to resolve this, we apply a 
nonuniform effective medium approximation that trans- 
forms networks with intermittent throats or significant 
variations in conductance into well-connected ones. 

The importance of connectedness within a pore net- 
work is highlighted by many results in the area of perco- 
lation theory . In particular, a very general result for 
such systems is that in the limit N — > oo, there is a crit- 
ical level of connectedness for the system which is called 
the "percolation threshold" - below this level, the overall 
network collapses into small, disjoint clusters, and hence 
there is no transport at the macroscopic level [27] . 

In general, it is extremely difficult to derive the per- 
colation threshold, or transport behavior in the vicin- 
ity of the threshold, for all but the simplest of network 
structures. This is especially pertinent to the proposed 
network model, since it has an irregular structure and 
correlations between the pores and throats. However, a 
simple analysis of the proposed model can identify some 
important structural properties. Using ty, the probabil- 
ity that a throat is present, we can derive the expected 
coordination number of a pore with extended coordinates 
y as [35] 

z(y) = (N - 1) / f{y )^>(y, y°)dy° (6) 
Js 

The value of z is crucial to overall transport in the 
medium - if this value is high, the network is well- 
connected and using the mean conductance for every 
bond should prove to be a reasonable approximation. Un- 
like with N, we cannot simply assume that z is large 
- pore networks extracted from 3D images of various 
porous media [lOj I28H32] have usually had average co- 
ordination numbers between 2 and 7. 

To illustrate the effects of low values of z, consider the 
simple case of a regular lattice pattern, in which every 
throat independently has probability p of being present 
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(.9 = .9o ) an d probability 1 — p of being absent (g = 0). 
Using a single-bond EMA, Kirkpatrick [15] has shown 
that the global transport properties of such a network are 
roughly the same as those of an identical lattice pattern 
in which every throat is present, with conductance [33J 



9 = 9oP 



z — 2p 



(7) 



As p — > 1 or z — > oo, Eq. ^ approaches the mean value 
approximation, g = gop. However, in cases where p and 
z are small, it is clear that the mean value approxima- 
tion is wrong by a significant margin, and hence a better 
approximation is required. 

Unfortunately, in general, single-bond EMA cannot ac- 
count for behavior near the percolation threshold of the 
system, because of a failure to account for the effects 
of correlations between the throats in the system [18] 
- in the terminology of percolation theory, the correla- 
tion length of the system diverges at the threshold [27] . 
Nonetheless, the accuracy of EMA tends to be quite good 
for systems that are not in the vicinity of their percola- 
tion thresholds [T5] . For the treatment of cases where z is 
small, with significant stochastic variation in the throats 
(especially when there is a high probability of absence) , 
we apply a generalization of single-bond EMA in which 
the bond conductances are not uniform, but instead vary 
depending on the extended coordinates of the pores they 
connect. The result of this approach is the transform- 
ing the probabilistic conductances gij into deterministic 
effective values g^. 

Unlike Kirkpatrick's approach for the regular lattice, 
we do not consider the effective conductance to be equal 
for every throat; instead, we account for the dependence 
of a throat's properties on the properties and locations of 
the pores it connects. Consequently, we solve for an an 
unknown function g^ = g(yi,yj) - effective throat con- 
ductance as a function of pore locations and attributes, 
creating an effective medium which is heterogeneous and 
anisotropic in nature. Note that the resulting effective 
medium is continuous in some sense - the pores are ef- 
fectively smeared out over their distributions, and the 
throat conductances are described by a function g. 

Although for regular lattices the effective medium can 
have low average coordination numbers, the randomness 
in pore properties and locations in the proposed model 
means that in our case the effective medium must usually 
have both large N and a large coordination number z, in 
order to properly eliminate the local fluctuations in the 
original network. In terms of the original network, the 
value z is the number of pores in what we call the locale 
of pore i - that is, the region of the medium to which pore 
i can (with nonzero probability) have direct connections. 
Fortunately, the structure of real world porous media is 
generally in line with this assumption - direct connections 
tend to occur only between nearby pores, but there are 
generally many other pores in the vicinity of any given 
pore. 

In order for the approximation proposed in this pa- 



per to be valid, the overall statistical properties of the 
medium (i.e. the distributions f(y) and f(g | y 1,1/2)) 
must not change significantly within a single locale. This 
type of assumption is typically almost inevitable in the 
derivation of a continuum model, although here it is only 
required on relatively small scales and so we can still ac- 
count for macroscopic heterogeneity such as spatial vari- 
ation in composition, structure or porosity. Aspects of 
the pores such as pore size can also be included in these 
statistical properties, but it is important to note that 
in such a case the locale consists of pores of varying 
sizes, and sharp changes in throat statistics between con- 
nected pores are allowed only with negligible frequency. 
In real-world porous media with fractal properties (coal 
may serve as a notable example) this type of approach 
should be entirely reasonable [HI 151] : in such cases, we 
expect transport to occur via a cascade process in which 
flow between large pores and much smaller ones occurs 
primarily via pores of intermediate size. 

The self-consistent condition of effective medium the- 
ory states that the fluctuations over a section of the orig- 
inal stochastic network, compared to its replacement in 
the effective medium, must average to zero. In our case, 
we consider this condition when applied to a single throat 
within the network connecting pores with extended co- 
ordinates taking the fixed values y, and y^; analogously 
to Kirkpatrick |15j . we express this condition as 



f{9 I Vi, Vj) — t, \ d 9 

9 + h(yi,yj) 



0, 



(8) 



where h(yi,yj) is the expected two-point conductance 
between pores at positions yi and yj , in the absence of 
a throat directly connecting them - this is the net flow 
through the network resulting from externally forcing 
a unit difference in concentration between those pores. 
This two-point conductance is usually defined in terms 
of Green's functions [35] : a detailed discussion with ex- 
act values for many infinite regular lattices can be seen 
in [35]. This equation allows the throat conductance to 
vary with the locations of the pores and the orientation 
and length of the throat, resulting in a heterogeneous and 
anisotropic effective medium. With this in mind, we term 
this approach the Nonuniform Effective Medium Approx- 
imation (NEMA). In its most general form, Eq. Q is 
very difficult to solve, since the values for h between any 
two pores depends on every other throat in the effective 
medium, and consequently on the values of the function 
g, while g in turn depends upon h. However, under cer- 
tain conditions, we can derive approximations for g and 
h which make the problem much more tractable. 

As long as almost all bonds in the original network have 
a very low probability of being present, $ < 1, which 
naturally results in z -C z and consequently (as needed) 
a well-connected effective medium, we can approximate 
9 by 



givuVj) = - 
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f{g I vuVj) 



1 9 + ^IKyuyj 



■ dg. 



(9) 
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For natural porous media, the assumption of $ < 1 
should prove to be a reasonable one, as there tend to 
be many pores in the vicinity of any given pore, but rel- 
atively few of these tend to be directly connected. For 
the details of the derivation of Eq. ([9]), see Appendix [B| 

It is clear that as long as an expression for the unknown 
function h(yi,yj) can be determined, Eq. ^ specifies 
the values of the conductances in the effective medium. 
Fortunately, the value of h is only important within a 
pore's locale, since long throats should be very unlikely 
and have very low conductances, resulting in negligible 
values of g from Eq. (Jf) ) regardless of h. Since the effective 
medium must have a large coordination number, we can 
see that in general h 3> g and so the presence or absence 
of a direct connection is insignificant for calculating h. 

To approximate the conductance between pore i and 
pore j, we consider it to consist of three conductances 
in series - two input conductances from the pores into 
their locales, which we represent by the unknown func- 
tion h ln (y), and a third term which comes from the 
interconnection between these locales. 

For neighboring pores, we expect that in a porous 
medium the locales of these pores should overlap signifi- 
cantly [37], resulting in 0(z 2 ) connections to one another, 
causing the third term to be negligible in comparison 
with the input conductances. This gives 
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(10) 



where we approximate h- m through connections of each 
pore to the distribution of other pores in the effective 
medium, i.e. 

kn(y) = (N - 1) / f(y°)g(y, v°) dy°, (n) 

Js 

whereas integration occurs over values of y° in the ex- 
tended space, with corresponding volume element dy° . 

By substituting Eq. ( 10 ) into Eq. ^ and the result 
into Eq. (Ill, and dividing both sides by h- m (y), we get 



(N-l)f(g,y°\y) 
sJq l + kn(y) [1/5 + VMi/°)] 



dgdy°. (12) 



This equation implicitly defines the function h- m (y), 
which can then be used in Eqs. (10 1 and Eq. ([9| to de- 



termine the values of g within the effective medium. 

Using the previous assumption that the distribution 
for conductances and properties does not change signifi- 
cantly within a given locale, and since h- ln primarily de- 
pends on the values of g, it follows that for y° within the 
locale of y, h in (y°) ps h in (y). With this approximation, 
the denominator of the integral becomes independent of 
y°; in conjunction with Fubini's theorem this gives 
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N- 1 



1 



hin(y)/g 



! f(g,y°\y)dy°dg. (13) 
Js 



The inner integral evaluates to the distribution of con- 
ductances coming out of a single pore with extended co- 
ordinates y, i.e. 



fig I y) = I f(9,y°\y)dy°, 
Js 



resulting in 



f(g I v) 



N-l J 2 + h in (y)/g 



dg. 



(14) 



(15) 



This gives the value of h ln for any y, which, together with 
Eqs. ( [To| and Eq. (J9|, specifies the throat conductances 
within the effective medium. 



IV. CONTINUUM APPROXIMATION 

A. Reduced master equation 

Applying the EMA described in the previous section 
accounts for the stochastic variation in throat conduc- 
tances in the original network. However, this still leaves 
a system with an equation for each individual pore, i.e. 



drrii 
~dT 



(16) 



which is merely the master equation Q but applied to 
the effective medium rather than to the original random 
network. In order to simplify this equation, we wish to 
transform it into a continuous form, wherein the amount 
of substance m and concentration c are considered over 
the extended space rather than within individual pores. 
This is achieved by considering the ensemble of realiza- 
tions of the pore positions and properties in the effective 
medium, and then converting to continuous form by tak- 
ing the limit of a summation over discrete regions. 

Assume that the extended space, S, is partitioned into 
countably many regions Sj (for J e N) such that each 
Si is a neighborhood of some point in extended space, 
yi. Here we mean "partition" in the usual mathematical 
sense, i.e. that the regions Si must be pairwise disjoint 
and cover S. If we let Aj be the set of pores contained 
within Si, i.e. Ai = {i | y^ € Si}, it follows that, in 
turn, the sets Ai form a partition of the set of all pores, 
A = {l,...,N}. 

The total amount of substance within the region Si is 
given by 



Mi=J2 



(17) 



and so, applying Eq. Q, the transport equation for Mj 
in the effective medium is given by 



dMi 
dt 



(18) 



() 



Since gij — g(yi, yj), and the Sj are neighborhoods in ex- 
tended space, it follows that as long as they can be made 
sufficiently small, with g continuous over each region (i.e. 
the set of discontinuities of g must have measure zero) , 
then for any i € Aj,j G Aj we have 



ho = g{yi,yj) 



(19) 



Since the value of the sum in Eq. ( 18 ) is finite, we can 
change the order of summation. If we denote the number 



of pores in the region Si by Ni, then Eq. (18) becomes 



dMi 
dt 



Jen V jeA., ieA z 



(20) 



As long as either (1) all pores have equivalent volume, 
Vb, or (2) the volume is included as a coordinate of 
the extended space (in which case we can set = In Vi 
without loss of generality) , the volume of each pore in a 
region can be considered equal, i.e. Vi = Vi. We now 
introduce averaged concentrations within each region, 



Ci = 



Mi 
~NiVi 



1 

]v7 



ieA 7 



(21) 



It is important to note that these concentrations are 
based on the total pore volume in each region, not the 
volume of physical space in each region. With these con- 
centrations, we can now write the master equation in 
terms of the regions, 



dt 



^giy^yMNjiCj-d). 



(22) 
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If we now consider the pores as being probabilistically 
spread over the space, the quantities Ni become the ex- 
pected number of pores in each region. These can be 
approximated using the marginal distribution f(y) as 



N T = N [ f(y)dy. 

J Si 



(23) 



If we also let the regions become infinitesimal, we get the 
reduced master equation, 



dM 
~di~ 



(y,t) = N 2 f g(y,y°)f(y)f(y°) 
Js 



[C(y°,t)-C(y,t)}dy°, (24) 



where C is the local concentration of the substance per 
unit pore volume, and M is the local concentration of 
the substance per unit of extended spatial volume. 



vanishes outside its locale, while the functions g, and / 
must not change significantly within this locale. Based 
on this, we expect C to behave similarly, and so we use 
power series to derive an approximation for Eq. ( 24 1 in 



the form of a diffusion-like Partial Differential Equation 
(PDE). 

We introduce Ay = y° — y = (Ay 1 , Ay 2 , . . .), and let 
gi satisfy 

s(y,y°) = .g 2 (y-y°,y + y°) = g 2 (Ay,2y + Ay). (25) 

Notably, since g is symmetric, it follows that gi must be 
an even function with respect to its first argument, Ay. 
The introduction of g 2 allows us to rewrite the reduced 
master equation as 

dM 



= N- 



[ g 2 (Ay, y + Ay)f{y)f(y + Ay) 
Js 



[C(y + Ay,t)-C(y,t)]d(Ay). (26) 



Taking power series expansions for 52 (Ay, y + Ay), 
f(y + Ay), and C(y + Ay,t) - C(y,t) yields a diffu- 
sional approximation - for detail on this derivation, see 
Appendix [C) Here we introduce 



(y) = Nf(y)V(y), 



(27) 



which means that M(y) = (f>(y)C(y). In the case where 
there are three spatial dimensions and no additional 
properties, i.e. y = (x 1 , x 2 , x 3 ), (j> is in fact the local 
(averaged) porosity of the medium. If <j) does not depend 
on time, then (using the Einstein summation convention 
over Greek indices) the final equation takes the form 



,dC_ 
~dt 



d 



(y 



\D a/3 ( 



y) 



dC 



(28) 



where the elements of the diffusion tensor D are given 
by 



y) 



Nf(y) 
W(y) 



Ay a AyPg 2 (Ay,2y)d(Ay). (29) 



It is important to note that while C is not conserved, 
Eq. (28 1 is conservative with respect to M = <pC, as 



should be the case for transport with no sources or sinks. 
The above integral is assumed to converge; any case 
where it does not would likely be a case of anomalous dif- 
fusion, which cannot be represented by the above equa- 
tions. It can be seen that this equation takes a form 
very similar to the Fokker-Planck (Kolmogorov forward) 
equation, which governs the PDF of an Ito-type random 
walk. However, the Kolmogorov forward equation dif- 
fers in that the diffusion coefficient is differentiated twice 
with respect to y. 



B. Diffusional approximation 



RESULTS 



In principle, the reduced master equation can be solved 
directly, but a simplified form is still desirable. Fortu- 



nately, we know from Section III that for any given y, g 



In this section, the accuracy of the approximation 
methodology described in Sections |III| and |IV| is evalu- 
ated for networks with different types of distributions for 
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conductances by comparison against numerical simula- 
tion results. For convenience, we abbreviate the over- 
all methodology as NEMA-DCA (Nonuniform Effective 
Medium Approximation and Diffusional Continuum Ap- 
proximation). As a further point of comparison, a more 
direct application of the diffusional approximation is 
given, in which each conductance is first replaced by its 
mean or expected value - this is precisely the simplified 
approach suggested at the very beginning of Section [Til] 
we term this approach the Averaged Diffusional Approxi- 
mation (ADA). Note that the use of the effective medium 
approximation by itself (i.e. without a continuum ap- 
proximation) does not directly yield useful results, since 
the resulting transport equations are not less complex 
than the original master equation. 

To simplify the mathematics, we consider cases where 
the conductances do not vary with space or direction, but 
have different probability distributions with respect to 
the throat length. Three characteristic types are consid- 
ered for the distribution of nonzero conductances within 
the medium: 



A single, discrete value (Section VD) 



A continuous distribution of conductances (Sec- 
tion V E I . For comparison, this distribution is cho- 



sen to have the same shape as in the previous case, 
but independently of throat length. 

• Conductances inversely proportional to throat 
length (Section [VF|. 



A. Simulation parameters 

To simplify the calculations, the pore networks are con- 
sidered to be contained within an n-dimensional cube 
having side length w, in which the pores are indepen- 
dently distributed with a uniform distribution. The num- 
ber of pores is set to N = 10000, and the pores are as- 
sumed to be identical in terms of their additional proper- 
ties £, resulting in all pores having the same volume, Vq. 
As such, the extended space can be considered to con- 
sist of only the spatial dimensions, y = x. Consequently, 
the pore network is contained within a hypercube of n- 
volume w n , and within this cube f(x) = w~ n , resulting 
in constant porosity 



y 



0o = NV /w n , 



(30) 



For further simplification, we consider a case where the 
probability distribution for the conductance of a throat is 
determined only by the quantity = | ] as* — ac_j ] | , which 
is the distance between the pores - this is directly related 
to the throat length, though in general they are not equal 
since the pores are not single points and the throats may 
be curved. Consequently, f(g \ Xi,Xj) — f(g | r^), 
and 52(Aa?, 2x) = g(r). As a result of this, the simulated 
networks can be described by a single diffusion coefficient 
which is independent of position and orientation - that 



is, they are homogeneous and isotropic with respect to 
macroscopic diffusivity. 

In order to localize the connections within the medium, 
we define the throat probability function as a rectangular 
function of the inter-pore distance, 



0, if r > r r - 
P, if r < r r , 



(31) 



From Eq. §6§, we can see that the expected coordi- 
nation number for any given pore is independent of its 
location, and takes the form 



P NK n 



pNK n r^ 



(32) 



where K n is the "surface area" of an n-sphere of unit 
radius. For convenience, we select as a characteristic 
distance tq the root mean squared distance between con- 
nected pores in the effective medium, which in this case 
is described by the equation 



(33) 



B. Analytical methods 

Given the simulation parameters above, application of 
the EMA described in Section IIIII results in a new net- 
work whose average coordination number is z = z/p, 
with no throats connecting pores further than r m apart, 
i.e. g(r) = for r > r m . On the other hand, for r < r m , 
there is a single value of h- m that is constant throughout 
the medium, and Eq. ( 15 ) takes the form 



1 



f(g) 



A - 1 Jo 2 + h in /g 



dg. 



(34) 



For r < r m , we have h = h ln /2 and so Eq. (|9|) takes the 
form 



9( r ) = 



f(g 



o 1/g + 2/hi 



dg, 



(35) 



while the diffusional approximation of Eq. ( 28 1 takes the 
simplified homogeneous and isotropic form 



dC 



x,t) = DV 2 C{x,t), 



in which the diffusion coefficient is 



D 



2K)C Jo 



g(r) dr. 



(36) 



(37) 



C. Numerical method 

In order to calculate numerical results for compari- 
son with the derived equations, a relatively simple ap- 
proach has been used. The graph is constructed as a 



symmetric N x N sparse matrix, in typical fashion for 
a weighted graph - throat conductances are represented 
by the nonzero values. In a modified form, this matrix 
can represent the master equation for every pore in the 
system as a matrix equation, which can be solved subject 
to appropriate boundary conditions. 

In these examples, steady-state solutions to the mas- 
ter equation are considered, in which the boundary con- 
ditions are specified by setting pores on opposing sides 
- specifically, within distance e of opposing faces of the 
cube - are fixed at different levels of concentration, and 
the resulting flow rate Q is calculated using the aforemen- 
tioned matrix method. This allows simulated diffusion 
coefficients to be calculated as 



}w(w — 2e) 
NV AC : 



(38) 



where AC is the difference in concentration between the 
two sides. 

The following sections will compare the numerically 
derived diffusion coefficients to analytical results using 
ADA (labeled as D a ) and using NEMA-DCA (D cm ). In 
order to display the nontrivial dependence of the diffusion 
coefficients on z, the diffusion coefficients shown in the 
figures are normalized to the dimensionless form 



D = D 



nV 
9o r l 



(39) 



where go is a characteristic conductance value. Although 
go is defined differently in some of the cases below, in each 
case the scaling factor nVo/go r o * s * ne same f° r both the 
analytical and numerical methods. 



D. Single- valued distribution 

In this class of networks is evaluated in which 

the conductance values, when present, always have a 
fixed conductance value, go- As stated previously, each 
throat has probability p — z/z of being present if the 
distance between the pores is no greater than r m , and 
zero probability of being present otherwise. 

Applying the ADA method would use g(r) = (z/z)go; 
implementing this instead of g in Eq. (37) predicts a dif- 
fusion coefficient of 



2 nVo ' 



(40) 



On the other hand, for the NEMA-DCA method, 



Eq. (34) gives 



9o(z-2). 



Applying Eq. ( 35 ) for r < r m leads to 

z-2 



-go, 



(41) 



(42) 




FIG. 1. (Color online) Results for the single-valued conduc- 
tance distribution. The circle, cross, and diamond markers 
are the results of simulations for n = 1, n = 2, and n = 3 
respectively, the dashed line is derived from ADA, and the 
solid line is derived from NEMA-DCA. 



which, when used in Eq. (37), gives 



2 nVo 



(43) 



A comparison of simulation results in 1, 2 and 3 di- 
mensions to the predicted diffusion coefficients for the 
single- valued conductance distribution is given in Fig. [T] 
Note that the axes, z and D, are both dimensionless. 
The graph demonstrates the accuracy of NEMA-DCA for 
varying numbers of physical dimensions, n, and varying 
values of z. While ADA predicts the asymptotic behav- 
ior as z — > oo, its relative error is rather high at lower 
values of z. On the other hand, NEMA-DCA performs 
quite well for z > 2. Nonetheless, neither approach cor- 
rectly predicts the behavior near the percolation limit, 
nor the percolation limit itself - this is to be expected, 
as EMA methodologies in general fail in approximating 
percolation. 

Although our networks take an irregular form, the re- 
sults shown in Fig. [I] are quite similar to those for Kirk- 
patrick's simple regular network (15| . The predictions of 
NEMA-DCA reflect the tendency for low coordination 
numbers to significantly reduce transport coefficients in 
physical porous media, as compared to the predictions of 
a simple diffusional approximation. 



E. Distance-independent inverse distribution 

As before, we consider n — 1, with no throats for pores 
further than r m apart, while throats covering shorter 
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FIG. 2. (Color online) Results for the distance-independent 
inverse conductance distribution, all for n = 1. The circle 
markers are simulation results, while the solid line is derived 
from NEMA-DCA. No line is shown for ADA, as it predicts 
infinite coefficients in this case. 



distances have a fixed probability z/z of being present. 
However, in this case, when a throat is present its con- 
ductance is determined (independently of r) by the dis- 
tribution 



this) 



-g 



(44) 



where k is a constant. As mentioned earlier, the star 
superscript indicates that this distribution is conditional 
on a throat actually being present, i.e. G > 0. For 
this example, we define the characteristic conductance 
as go = k/r . Notably, ADA would fail entirely in this 
case, as g = oo and so 



(45) 



However, the application of NEMA removes this sin- 
gularity. By using Eq. (44 1 in Eq. ( 34 ) , we get 



kz 



1 



dg. 



Ik/r m 2g 2 + h n g 
If we define the function B implicitly by 

'2B(s 



1 + B(s) = exp 
then hi n can be expressed as hi n = 2kB(z)/r ri 



(46) 



(47) 



However, in this case Eq. (|35|) becomes 



2kB(z) 



(48) 



FIG. 3. (Color online) Results for the inverse-distance con- 
ductance distribution, all for n = 1. The circle markers are 
simulation results, the dashed line is derived from ADA, and 
the solid line is derived from NEMA-DCA. 



for r < r m . Note that, as in the single- valued case, this 
value does not depend on r, and so every throat in the 
effective medium has the same conductance. 



Using Eq. ( 48 ) in Eq. ( 37 1 we derive the diffusion co- 
efficient 



B{z) g rl 
V ■ 



(49) 



A comparison of simulation results to the predicted 
diffusion coefficients for the distance-independent inverse 
conductance distribution is given in Fig. [2j Once again, 
NEMA-DCA performs quite well, though it does not pre- 
dict the percolation threshold accurately. On the other 
hand, ADA in this case is only valid for infinite z, as it 
predicts infinite diffusion coefficients. 

Although this example presents an extreme case, it 
serves to illustrate the use of effective medium techniques 
in order to provide a realistic account for the effects 
of sparse high-throughput connections within irregular 
porous media. Such connections can be present in physi- 
cal porous media, but it is clear that the overall through- 
put would nonetheless be bottlenecked by the smaller 
connections, as is predicted by NEMA-DCA for lower 
values of z. 



F. Inverse-distance distribution 

For this example, we once again have throats being 
present with probability zjz'xir < r TO , and as in the 
previous example we select n = 1. However, in this case, 
whenever such a throat is indeed present, its conductance 
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takes the value g = k/r. For simplicity and ease of com- 
parison to the previous case, we define the characteristic 
conductance as g$ = k/r Q . 

In this case, the ADA approach would use g(r) — 
(z/z)(k/r), which when used in Eq. (37) gives 



D n = 



V3 

~4 



V 



(50) 



To apply NEMA-DCA, we require a PDF for the con- 
ductances. Since f*(r) — l/r TO , it follows that 



/*(<?) 



9 



(51) 



Notably, this is, in fact, the same distribution of con- 
ductances as in the previous example, and the inter-pore 
distances r also retain the same distribution as before, 
but this case is quite different from the previous one due 
to the correlation between the conductance and inter- 
pore distance. Consequently, if the dependence of con- 
ductance on distance were ignored, NEMA-DCA would 
make the same prediction here as in the previous case. 

Since the conductance distribution is the same, Eq. [34] 
takes an identical form to that in the previous case, and 
so once again we have h m — 2kB (z) /r m , while the func- 



tion B(z) is determined by Eq. g7j). With Eq. ([35]) for 
r < r m this gives 



g( r ) = 



z/z 



r/k + r m /[kB(z)Y 



(52) 



and, from Eq. (37), 



kz 
2Vor, 

4 



z-2 



B(z) 



■rn/B(z) 

9o r o 

Vn 



dr. 



(53) 



A comparison of simulation results to the predicted 
diffusion coefficients for the inverse-distance conductance 
distribution is given in Fig. [3] As can be seen, NEMA- 
DCA is very accurate except near the percolation limit, 
while ADA predicts the asymptotic behavior and is accu- 
rate only for larger z. Like the first example, this example 
illustrates the need to account for the effects of limited 
interconnectedness in order to achieve a realistic model. 



G. Overall discussion 

In contrasting the last two cases, with differing inverse 
distributions, it is clear from Figs. [2] and [3] that the re- 
sults are very different - the radius-independent distri- 
bution gives a higher diffusion coefficient by a factor of 
2 at z = 10, and it is clear from the shape of the curve 
that this divergence persists for increasing z. This shows 
the importance of properly representing various property 
correlations, especially between pores and throats, in the 



porous medium. Such correlations are the only difference 
between the cases in Fig. [2] and Fig. [3] and it is clear that 
these kinds of correlations are an essential aspect of real- 
world porous media. Unlike NEMA, a more typical EMA 
(e.g. Kirkpatrick [TS]) would not be able to account for 
this type of behavior. 

However, it is interesting to note that in each case, 
the NEMA predicts a percolation limit at z = 2, just 
as in Kirkpatrick's paper [T5] - this seems to be a uni- 
versal property of single-bond effective medium approx- 
imations, regardless of whether the network in question 
is regular or irregular in nature. Notably, as is demon- 
strated by the second example, the behavior away from 
z = 2 does not need to be linear. Although the numeri- 
cal results show that the behavior very close to the true 
percolation limit is not correctly predicted, as is gener- 
ally the case for effective medium approximations, it can 
be seen that the NEMA-DCA prediction starts to match 
the numerical results very closely at a short distance from 
the percolation limit. 

Moreover, the general interdisciplinary tool developed 
in this paper can be applied to Fickian diffusion, Knud- 
sen diffusion, and incompressible Newtonian Stokes flow; 
example throat conductance expressions for some simple 
throat geometries are presented in Table 1 of Appendix 

m 



VI. CONCLUSIONS 

In this paper, we have presented a general analytical 
methodology for analyzing various transport phenomena 
in irregular porous media. We have developed a gener- 
alized effective medium approximation which allows the 
conductances of the throats in the medium to vary de- 
pending on the properties and positions of the pores they 
connect. This is used in combination with a continuum 
approximation in the form of a diffusion-like equation to 
derive an explicit analytical expression for the transport 
coefficients, allowing for an efficient model of transport 
that retains the ability to account for microscopic effects. 

These analytical expressions have been tested against 
several numerical simulations with different types of dis- 
tributions for throat properties, demonstrating very good 
accuracy for the proposed methodology, except in the 
vicinity of the percolation limit. Furthermore, the re- 
sults demonstrate the ability of this method to account 
for several physical effects, including reduced transport 
due to low coordination numbers, and the effects of cor- 
relations between local properties of the medium. 
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Appendix A: Expressions for throat conductances 



Recalling that the master equation 

N N 



drrii 
~dT 



3=1 3=1 



(Al) 



applies to a variety of different types of transport, Table|T] 
gives examples of throat conductances for a variety of 
different cases. 

Case Expression for g Ref. 

Small orifices Long, narrow throats 





(L = 0) 


(L > a) 




A 


2Da 


Dna 2 /L 


[38] 


B 


vita 2 /A 


2vna 3 /(3L) 


m 


C 


a 3 /(3 M ) 


7ra 4 /(8/iL) 


m 


D 




ana 2 /L 





TABLE I. Expressions for throat conductance g for various 
transport phenomena. Where applicable, a is the radius of 
the throat, L is the length of the throat, D is the diffusion 
coefficient, v is the mean molecular speed, fi is the fluid vis- 
cosity, and a is the electrical conductivity of the material. 

The cases considered in Table U are: 

(A) Fickian Diffusion. 

(B) Knudsen Diffusion. 

(C) Incompressible Newtonian Stokes flow. Since 
the flow is steady, there is no time derivative term; Cj 
is interpreted as the pressure in each pore, while qij 
is the volumetric flow rate through the throat from 
pore i to pore j. 

(D) Flow of electrical current. Once again, the flow 
is steady; Ci is the voltage at node (pore) i, and qtj is 
the current flowing through the branch (throat) from 



node i to node j. Note that in the case of electrical 
transport, the "pores" and "throats" usually consist 
of the solid matrix rather than the void space. 



Appendix B: Approximation for g 

Starting with the effective medium criterion [Eq. 



f{9 I Vi, Vj) f7 dg = 0, 

9 + h{yi,yj) 



we apply Eq. ([2]), resulting in 



(l-*)f = ¥/ r{g\y i ,y 3 ) 9 -—ldg, 
h Jo g + h 



(Bl) 



(B2) 



where Vt, g and h are, as before, functions of yi and yj - 
the arguments have been dropped for compactness. This 
can be rearranged to give 



* r r($\vuvi) 



where 



x = i-* 



X Jo l/g + 1/h 

f*(g I yt,yj) 
o i + g/h 



dg, 



1 - 



dg 



(B3) 



(B4) 



Since /* is a PDF and the conductances must take pos- 
itive values, the integral within the brackets is bounded 
below by and above by 1, and so 1 — ^ < x < !■ For 
f < 1, we get x ~ 1 an d hence Eq. (B3) can be written 

as 



f(g I yuvj) 

l/g + l/h{y l ,y j 



dg, (B5) 



since the 8(g) component of the integral vanishes due to 
the l/g term in the denominator. This completes the 
derivation of Eq. Q . 



Appendix C: Derivation of the diffusional approximation 



Based on the assumption in Section 111 that the network has no long-range connectivity and its statistical properties 
do not change significantly within any locale, we expect that power series expansions should offer reasonable approxi- 
mations for / and g~2, and, based on the form of Eq. (26 1, for C . The power series expansions for C(y + Ay, t) — C(y, t), 
f(y + Ay), and g 2 (Ay, y + Ay), are as follows: 



8C 1 f) 2 C 

C(y + Ay, t) C(y, t) = Ay^ — (y, t) + -Ay" Ay* 9 (y , t) + 



(CI) 



f(y + Ay)=f(y) + Ay a -^(y) + .. 



(C2) 
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and 

g 2 (Ay, 2y + Ay) = g 2 (Ay, 2y) + ±Ay a ^(Ay, 2y) + . . . . (C3) 

Using these expansions while discarding third- and fourth-order terms allows the reduced master equation to be 
expressed as 

9M iV 2 f d 2 C f2 d~g 2 dC , n£ df dC 



ot - 2 j s ^^{^a^ + ^WW + 2f ^ h W d(Ay) ' (C4) 



where the arguments for the functions are dropped for 
compactness; as before, / and its derivatives are eval- 
uated at y, g 2 and its derivatives are evaluated at 
(Ay,2y), and C and its derivatives are evaluated at 

(y,t). 

Note that the first-order (drift) terms in this equation 
are zero. This follows from the fact that f(y) and C(y, t) 
are independent of Ay while g 2 is even with respect to 
Ay, resulting in 



Ay a g 2 d(Ay)=0, 



(C5) 



Using the product rule, as well as the aforementioned 
independence, Eq. (C4) can be expressed as 



DM _ d 
Ik ~ dy* 



AT 2 

IT 



f [ Ay a Ayf ) ~g 2 d(Ay) 
Js 



dC 



(C6) 



completing the derivation of the diffusional approxima- 
tion. 
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